#ifndef RATBEZIER_H
#define RATBEZIER_H

#include <cassert>
#include <algorithm>
#include <vector>

#include "Curves/Interfaces/icurve.h"
#include "Curves/CurveUtilities.h"

template< class IPoint >
class RatBezier : public ICurve< IPoint >
{
protected:
    /* Fields. */
    int degree;
    vector< double > weights;

public:
    /* Constructors. */
    RatBezier( );
    /* Destructors. */
    ~RatBezier( );

    /* Methods. */
    void AddPoint( const IPoint& point, double weight );
    void AddPoint( const IPoint& point, int point_index, double weight );

    void RemovePoint( );
    void RemovePoint( int point_index );

    IPoint Eval( double t ) const;
    vector< IPoint > ComputePoints( int n );

};

template< class IPoint >
RatBezier< IPoint >::RatBezier( )
{
    degree = -1;
}

template< class IPoint >
RatBezier< IPoint >::~RatBezier( )
{
}

template< class IPoint >
void RatBezier< IPoint >::AddPoint( const IPoint &point, double weight )
{
    ICurve< IPoint >::AddPoint( point );
    weights.push_back( weight );
    degree++;
}

template< class IPoint >
void RatBezier< IPoint >::AddPoint( const IPoint &point, int point_index, double weight )
{
    ICurve< IPoint >::AddPoint( point, point_index );

    vector< double >::iterator it = weights.begin() + point_index;
    weights.insert( it, 1, weight );

    degree++;
}

template< class IPoint >
void RatBezier< IPoint >::RemovePoint( )
{
    ICurve< IPoint >::RemovePoint( );
    degree--;
}

template< class IPoint >
void RatBezier< IPoint >::RemovePoint( int point_index )
{
    ICurve< IPoint >::RemovePoint( point_index );

    vector< double >::iterator it = weights.begin( ) + point_index;
    copy( it + 1, weights.end( ), // Source
          it );                   // Destination
    weights.pop_back( );

    degree--;
}

template< class IPoint >
IPoint RatBezier< IPoint >::Eval( double t ) const
{
    assert( 0.0 <= t && t <= 1.0 );

    IPoint P[degree + 1];
    double W[degree + 1];
    double s;

    for ( int i = 0; i <= degree; i++ )
    {
        P[i] = this -> control_points[i];
        W[i] = weights[i];
    }

    for ( int k = 1; k <= degree; k++ )
        for ( int i = 0; i <= degree - k; i++ )
        {
            W[i] = ( 1.0 - t ) * W[i] + t * W[i+1];
            s    = t * W[i+1] / W[i];
            P[i] = ( 1.0 - s ) * P[i] + s * P[i + 1];
        }

    return P[0];
}

template< class IPoint >
vector<IPoint> RatBezier< IPoint >::ComputePoints( int n )
{
    if ( this -> computed_points_up_to_date &&
         this -> computed_points.size( ) == static_cast< unsigned int >( n ) )
        return this -> computed_points;

    this -> computed_points.resize( n );

    for ( int i = 0; i < n; i++ )
        this -> computed_points[i] = Eval( static_cast<double>( i ) / ( n - 1 ) );

    this -> computed_points_up_to_date = true;

    return this -> computed_points;
}



#endif // RATBEZIER_H
